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ABSTRACT 



In this thesis, an adaptive two dimensional least mean squares (LMS) algorithm and 
a recursive least squares (RLS) algorithm are developed from the one dimensional algo- 
rithms. Design of the two dimensional LMS and RLS algorithms are studied for accu- 
racy based on the results of a two dimensional system identification model which was 
used in testing the algorithms. Application of the two algorithms is demonstrated 
through computer simulation in which the adaptive filters are employed in a noise 
canceler and an adaptive line enhancer and applied to an image processing problem. 
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I. INTRODUCTION 



The area of digital signal processing has experienced a rapid growth in the last dec- 
ade. Reasons for this have been the tremendous advances in integrated circuit technol- 
ogy, and some significant developments in digital processing techniques achieved during 
this period. Included in these developments are methods which extend certain one di- 
mensional digital signal processing techniques to two dimensions. This extension is by 
no means trivial. Three significant factors must be considered: 

1. more degrees of freedom are available in a two dimensional system; this gives a 
system designer more flexibility than the one dimensional case, 

2. one dimensional problems generally involve considerably less data than two di- 
mensional problems, and 

3. the mathematical methods for handling two dimensional systems are generally less 
complete than those for one dimension. 

As the techniques for processing multidimensional data have improved, the appli- 
cations of digital signal processing have spread from one dimensional to two dimen- 
sional, to n-dimensional data. There exists many physical phenomena that inherently 
depend on two or more independent variables. In the prediction of weather and in 
seismic analysis, the data depends on more than one independent variable. Moreover, 
data originating from one dimensional processes can, in some cases, be considered two 
dimensional. For instance, data from periodic or cyclic processes can be represented as 
two dimensional arrays by using their periodicity. 

Besides the general applicability of two dimensional signal processing in the above 
cases, other areas which have experienced significant growth in recent years include ra- 
dar, sonar, and radio astronomy. The two dimensional processing of images is also a 
very important one. Images depend on two spatial variables, and are continuous in 
nature. However, if we digitize them and assume linear models in their formation, dis- 
tortion, and recording we then have techniques that can be used in their processing. 
Depending on the applications, different processing techniques are used, notably: en- 
hancement and restoration of images, and segmentation and encoding of images. For 
details, see References 1,2,3, 4 . 

This brief discussion about applications of two dimensional digital signal processing 
shows that they are found in a wide variety of fields. In this thesis, we are interested in 
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extending one dimensional adaptive filtering techniques to two dimensions and then 
applying this in the area of image restoration. Once the transition from one to two di- 
mensions is understood, the extension to n-dimensional signal processing is fairly 
straightforward. 

A. OBJECTIVES OF THE THESIS 

The use of the Wiener filter has proven to be a very powerful tool in the area of 
image restoration [Refs. 1,2, 4 ]. However, one disadvantage is the fact that the Wiener 
filter operates under the assumption that the image is stationary which generally is not 
the case in an image processing problem. In one dimensional signal processing the 
Wiener filter concept has provided a basis for various adaptive filtering algorithms. 
Within these adaptive algorithms, the filter possesses characteristics which can be mod- 
ified to achieve some end or objective and is usually assumed to accomplish this "adap- 
tation" automatically without the need for substantial intervention by the user. The 
adaptive filter can "learn" the signal characteristics when first turned on and thereafter 
can track changes in these characteristics. The first objective of this thesis is to examine 
a two dimensional Wiener model for image restoration which can then be extended to 
adaptive algorithms. Due to its relatively low computational requirements and the fact 
that it will work in a variety of signal environments, the least mean square (LMS) algo- 
rithm will be investigated first. 

The second objective of this thesis is to develop a second two dimensional adaptive 
algorithm. In this case, we will work with the recursive least square (RLS) algorithm 
which offers faster convergence than the gradient-search-type algorithms but generally 
involves a greater cost per data sample and more numerical difficulties. 

The final objective is to implement the LMS and RLS algorithms within a system 
identification model, a noise canceler, and an adaptive line enhancer. Each algorithm 
possesses several variants and various parameters which can be modified. A represen- 
tative sample of the various outputs will be examined and the results compared in order 
to see w'hich provides the optimum solution under given conditions. 

B. ORGANIZATION OF THE THESIS 

Chapter II is designed to review the development of the one dimensional Wiener 
filter and then extend it to two dimensions where it can be incorporated into a two di- 
mensional LMS adaptive algorithm. Computer simulation of the algorithm within a 
system identification model is performed and the results are shown. The two dimen- 
sional RLS algorithm is derived in Chapter III and again the results of the algorithm 
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within a computer simulated system identification model are shown. Chapter IV con- 
tains the results of implementing the LMS and RLS algorithms in a noise canceler and 
an adaptive line enhancer. Conclusions concerning the results are also presented. 
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II. TWO DIMENSIONAL ADAPTIVE LEAST MEAN SQUARE 

ALGORITHM 



In this chapter, we will review the derivation of the one dimensional Wiener filter 
and then develop an algorithm to extend it to the two dimensional case. We will use the 
two dimensional Wiener filter to achieve a two dimensional least mean square (LMS) 
adaptive filter algorithm which will be demonstrated to be useful in image processing by 
applying the algorithm in a noise canceler mode and in an adaptive line enhancer con- 
figuration for the restoration of images corrupted by noise. 

A. ONE DIMENSIONAL WIENER FILTER 

The problem of estimating one signal from another is one of the most important in 
signal processing. In many applications, the desired signal is not available or observable 
directly. Instead, the observable signal is a degraded or distorted version of the original 
signal. The signal estimation problem is to recover, in the best way possible, the desired 
signal from its degraded replica. One typical example which we will deal with in this 
paper is an image recorded by an imaging system that has been corrupted by noise. The 
problem is to undo the noise induced distortion and restore the original image. 

This represents the classic one dimensional problem in communication theory where 
we must obtain an estimate of a signal of interest, which can be observed in the presence 
of some additive noise. In other words, the available information about the signal, 
s(n), is contained in the received signal: 

u(n) = s(rt) + vv(rt) (2.1) 

where w(«) is the noise. Therefore, we must process this available signal u(n) through 
an optimal processor that produces the best possible estimate of s(w). 

In order to establish a two dimensional Wiener filter, we must first develop a one 
dimensional algorithm. This task has been approached from many different directions 
[Refs. 5,6,7]. The formulation by Haykin [Ref. 8] provides the most logical extension to 
two dimensions. First, we consider a tapped delay line filter similar to Figure 1 on page 
5. The filter consists of a set of delay elements, a corresponding set of adjustable tap 
gains or coefficients /i(l), h(2),..., h(M) connected to the tap inputs, and a set of adders 



4 




Figure 1. Tapped Delay Line Filter 

For summing the resultant outputs. The filter is driven by a random time series produc- 
ing the sequence //(«), u(n — 1),..., u(n — M + 1) as the M tap inputs of the Filter. 

We can express the signal produced at the Filter output, y(n), by the convolution 
sum: 



>•(") 




h(k)u{n-k+ 1). 



( 2 . 2 ) 



We desire a Filter which in some way minimizes the difference between some desired re- 
sponse, d(n), and the corresponding value of the actual Filter output. This difference can 
be denoted as 



e(n) = d(n ) - y(„) (2.3) 

where e(n) is called the error signal. In Wiener theory, the Filter is optimized by mini- 
mizing the mean-square value of this error signal, e(n) . 
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Let the mean-square value of the error be denoted by 

MSE=E{e\n)} (2.4) 

where E {.} is the expectation operator. This mean-square value is a real and positive 
scalar quantity, representing the average normalized power of the error signal, e(n). 
Substituting Equation (2.3) into (2.4) yields 

MSE = E {/(«)} -2 E M«Mn)) + E ly 2 («)}- (2.5) 

Next, substituting Equation (2.2) into Equation (2.5) and then interchanging the orders 
of summation and expectation in the last two terms, we get 



M 

MSE = E {d\n)} - 2 V h{k) E { d{n)u{n - k + 1)} 

M M 

+ ^ ^ h(k)h(m ) E {u(n — k + 1)u(k — m + 1)} . 

k=\ m = 1 



( 2 . 6 ) 



Assuming that the input signal u(n) and the desired response d(n) are jointly stationary, 
the three terms on the right-hand side of the above equation may be interpreted as fol- 
lows: 



1. The expectation E {^(w)} is equal to the mean square value of the desired response 
d(n): 



P d =E{d\n)}. 



(2.7) 



2. The expectation E {d{n)u(n — k -I- 1)} is equal to the cross-correlation function of 
the desired response d(n) and the input signal u(n) for the lag of k-1. We can 
therefore write the single summation term on the righthand side of Equation (2.6) 
as follows: 




h{k)E{d{n)u{n -*+!)} = 




h(k)p(k - 1) . 



( 2 . 8 ) 



3. Finally, the expectation E {u(n — k + l)u(n — m + \)} is equal to the 
autocorrelation function of the input signal u(n) for the lag of m-k: 



r(m — k)= E {u(n — k + l)w(« — m -1- 1)} . (2.9) 



Accordingly, we can rewrite the double summation term on the righthand side of 
Equation (2.6) in the form 




M 

y h(k)h(m) E{u(n — k + 1 )u(n — m + 1)} 

OT=rl M M 




y/ 2 {k)h{m)r{m - k ) . 

m=\ 



( 2 . 10 ) 
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Thus, substituting Equations (2.7), (2.8), and (2.10) into Equation (2.6), we find that the 
expression for the mean square error may by rewritten in the form 

M MM 

MSE= P d -2^h{k)p{k-\)+ '^h(k)h(m)r(m — k ) . (2.11) 

By differentiating Equation (2.11) with respect to h(k) and setting it equal to zero, we 
have the following set of M simultaneous equations 



M 

p{k - 1) = ^Ao(w)r(m - k), /c = 






( 2 . 12 ) 



This system of M simultaneous equations is called the normal equations with optimum 
filter coefficients as the unknowns. With the following definitions, 



W) 

W) 



Ao(A/) 



P( 0) 
/>(!) 



p(M-\) 



r(0) r(l) 

r(l) r(0) 



r{M- 1) 
r{M - 2) 



R = 



r(M — 1) r(M — 2) r(0) 



(2.13) 



(2.14) 



(2.15) 
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we can rewrite the normal equations of Equation (2.12) in matrix form as 



p = R h 0 . 



( 2 . 16 ) 



This represents the discrete-time version of the Wiener-Hopf equation. 

B. TWO DIMENSIONAL WIENER FILTER 

The following derivation parallels the development by Hadhoud and Thomas, how- 
ever this research and simulation were completed separately and prior to the publication 
of Reference 9. In order to be applicable for an image processing problem, we must 
extend the formulation in the previous section to two dimensions [Refs. 10,11]. This is 
accomplished by developing a basic two dimensional Wiener filter as shown in 
Figure 2 on page 9. Within this filter, we use two input images, the reference array X 
and the primary input array D. The primary input image D is a two dimensional array 
which represents the ideal image plus additive noise, while the reference image X is noise 
that is assumed to be correlated to the noise in the primary input. Both the input arrays 
are MxM in dimension. The Wiener filter is an NxN causal FIR filter with a set of 
weights Wj defined as 



which minimize the mean of the squared error, e p between the filter output and the de- 
sired input D. We designate j as the iteration number given by j = mM 4- n where m and 
n take on the values from 0 to M-l. Just as in the one dimensional case, the filter out- 
put, y(m,n), is the convolution sum of the filter mask and the reference input X which is 
given by 



Wj( 0,0) Wj( 0,1) ... Wj(0,N — 1) 

Wj, 1,0) Wj( 1,1) ... 1) 




(2.17) 



Wj(N — 1,0) Wj(N — 1,1) ... Wj(N — \,N — 1) 



v-i v-i 



y(m,n) = Z Z W A‘' k ) X ( m ~ /,« - *)• 



( 2 . 18 ) 



/=o E=o 



:0 



During the jth iteration the input from array X is represented by X t where 
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Figure 2. Two Dimensional Wiener Filter 
X(m,n) X(m,n— 1) 



A} = 



X(m,n — N + 1) 



(2.19) 



X(m -N+l ,n) X(m -N+\,n-\) ... X{m -N+\,n-N+\) 

/ 

Using the Equations (2.17) and (2.19), Equation (2.18) can now be written as 



9 



( 2 . 20 ) 






A'-l A'-l 

y Y.ivjiU) Xjm) 

% 0 k=0 



for the jth iteration. 

This output y{m,n ) can now be subtracted from one pixel D(m,n) of the array D to 
produce the error signal at the jth iteration 



A'-l A'-l 



ej = D(m,n) - > > Wj(l,k)X(m - l,n - k). 



( 2 . 21 ) 



s=o £=o 



Since the purpose of the Wiener filter is to provide a set of weights which minimize 
the MSE, w'e can denote it as 



MSE=E{ej) 



( 2 . 22 ) 



where E{.} is the expectation operator. 

Using a mathematical derivation similar to the one dimensional case of the previous 
section, we can see that 



A'-l A'-l 



e] = D\m,n) -2 Y Y \V{l,k) D(m,n ) X(m - !,n - k) 

/V— 1 am l » » 

+ Z Z Z X *(»> -l,n-k)X(m - p,n -q). 



(2.23) 



/=0 / c =0 p =0 <?=0 



Substituting Equation (2.23) into Equation (2.22) yields 



A— 1 A r -1 



MSE = E[Z) 2 (m,«)] - 2 > ) Wj{l,k) E\_D{m,n) X{m — l,n — A)] 



A'-l A'-l A'-l A'-l 

+ E E E E - i ’ n - *) - P'" - • 

/=0 £=0 p = 0 <7=0 



7=0 £=0 



(2.24) 



Defining P as the crosscorrelation matrix between the desired response D(m,n) and the 
reference input, R as the input autocorrelation matrix, and W as the optimum Wiener 
weight matrix, we can rewrite Equation (2.24) as 
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(2.25) 



A r -1 A’— 1 



MSE = E[D\m,n )] - 2 ^ Y JF/U) P(/,£) 

A— 1 A’-l A— 1 A— 1 /=0 

IIII Wj{l,k) Wj(p,q) R(p - l,q - k) . 



1=0 k=Q p=0 q = 0 



Finally, if we minimize the MSE with respect to WfJ,k) then we have 



N—\ N-\ 



P(W ' V jM R (P ~ U “ *)• 



(2.26) 



p = 0 <y=0 



This equation is the two dimensional equivalent to Equation (2.12). The matrix form 
of Equation (2.26) is given as 



P = RW 

where the elements of this equation are defined as follows 



w ra [/y - 

c*i] VU M - 



IR — AM-1 ] IR — A r +2^ IRol 



^0 



w = 



w 



N—\ 



P 0 
P\ 



(2.27) 



(2.28) 



(2.29) 



(2.30) 



P 



A'-l 
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Within the R matrix in equation (2.28) each element is a block Toeplitz matrix repres- 
ented by the equation 

R{p — l,q — k) = E\_X{m — l,n — k)X{m — p,n — q)~\ (2.3 1) 

where l,k,p,q range from 0 to N-l. 

C. TWO DIMENSIONAL LEAST MEAN SQUARE 

One means of obtaining an approximate solution for the optimum weight matrix, 
W, is the use of a two dimensional LMS algorithm which is depicted in Figure 3 on page 
13. This differs from Figure 2 on page 9 in that the error is used to update the filter 
coefficients before shifting the data window X } across the reference input for the next 
iteration. As in the one dimensional LMS algorithm, we are updating the coefficient 
matrix by adding the present matrix to a change proportional to the negative gradient 
of the error where the one dimensional instantaneous estimates of the gradient vector 
are based on sample values of the input and the error signal e{n)u(n) [Ref. 8]. For the 
two dimensional jth iteration, we define the updated matrix as 

Wj+^lVj-pejXj (2.32) 



where 

Wj+t updated weight matrix 
W } previous weight matrix 

p scaler multiplier controlling the rate of convergence and filter stability 
{ej)(Xj) estimate for the 2-D instantaneous gradient 

The previous equation can also be written as 



W J+ ,(l,k) = Wj(U ] ) + 2 p ( ej ) X(m - l,n - k). (2.33) 

These two equations give the two dimensional weight updating algorithm for the LMS 
adaptive filter. The algorithm we have developed here may be implemented without any 
form of matrix operations, averaging, or differentiation. 

The value of p may be chosen based on the desired tracking ability, steady-state 
mean square error, and convergence speed. In signal processing, there are several 
methods for determining a suitable value, however in many of these cases it requires 
knowledge of the eigenvalues and eigenvectors of the autocorrleation matrix. In image 
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Figure 3. Two Dimensional Least Mean Square 

processing, one method which does not require these values for the computation of n is 
trial and error based on output image. 

An alternate method used in one dimensional design which again does not require 
a priori knowledge of the autocorrelation matrix, discussed by Bitmead and Anderson 
[Ref. 12], is the normalized LMS. Using our previous notation, this method can be ex- 
tended to two dimensions by first considering the two dimensional LMS update equation 
(2.32). In this equation, we redefine the step-size parameter, n, for a given / and k as 
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nm = 



a 



(2.34) 



P + o](l,k) 



in which a, represents the normalized step size chosen between zero and two, /? is an- 
other small positive constant, and a 2 (l,k) is one of the values from the matrix 



»>, 0) oj(0,l) ... o](0,N-l) 

Oy(l,0) Oj (1,1) ... 



aj(N — 1,0) <Tj(/V— 1,1) ... o]{N- \,N- 1) 



(2.35) 



The element required from the matrix depends upon the current values of / and k being 
used for equation (2.33). The matrix may be be intitialized with the values of 
{X(m — l,n — k)) 2 and to update the values in o 2 , we use the following equation 

o] +l (l,k) = p(°](l,k)) + (1 - p)(X(m - l,n - k)f (2.36) 



where p is a weighting factor between zero and one. This ensures that the value of p 
does not become large enough to cause the algorithm to become unstable. 

D. IMPLEMENTATION 
1. System Identification 

In order to initially test the two dimensional least mean square algorithm, we 
established a system identification model shown in Figure 4 on page 15. Within this 
model the output of the known FIR Filter, d{m,n), is defined as 



d(m,ri) = .4 W(m,n) + .6 W(m — 1,/?) 
— .3 W(m,n — 1) + .2 W(m — 1 ,n — 1) . 



(2.37) 



The output of the two dimensional least mean square filter, y(m,n) , consists of a set of 
adjustable coefficients and is defined as 

y(m,n) = AO W(m,n) + A 1 W(m — 1,/?) 

+ A2 W(m,n — 1) + A3 W(m- \,n- 1) . (2 ' 38) 



The adaptive filter output y(m,n) is then compared with the known system output 
d{m,ri) to produce an error signal e{m,n), defined as the difference between them. 
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error 



Figure 4. Two Dimensional System Identification 

The operation of the adaptive filter is to minimize the error signal e(m,ri) by 
providing an adaptive process for adjusting the coefficients. For this adaptive process 
we use the update equation (2.32) for the two dimensional least mean square algorithm 
developed in the previous section 



IV J+1 = IVj-nejXj (2.32) 

For this model a 32x32 white gaussian noise matrix was used as the input and 
the rate of convergence can be be seen in Figure 5 on page 16 and Figure 6 on page 
17. Following 600 iterations all the coefficients had converged to within 10~ 3 of the ac- 
tual coefficient value and the error was 3;d0 -s . A computer program for this system 
identification model is given in Appendix A. 

2. Adaptive Noise Canceler 

The usual method of estimating a signal corrupted by additive noise is to pass 
the composite signal through a filter that tends to suppress the noise while leaving the 
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Figure 5. LMS System Identification-Rate of Convergence 

signal relatively unchanged. The noise canceler and adaptive line enhancer developed 
by Widrow [Ref. 13] are well-documented ways of doing that. Adaptive noise canceling 
is a variation of optimal filtering that is highly advantageous in many applications. It 
uses an auxiliary or reference input derived from one or more sensors located at points 
in the noise field where the signal is weak or undetectable. This input is filtered and 
subtracted from a primary input containing both signal and noise. The reference input 
and the noise in the primary input are therefore correlated, and as a result the primary 
noise is attenuated or eliminated by cancellation. 
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Figure 6. LMS System Identification-Rate of Convergence 

The basic noise canceling concept is illustrated in Figure 7 on page 18. A signal 
is transmitted over a channel to a sensor that receives the signal plus noise, n 0 . The 
combined signal and noise s + n 0 forms the "primary' input" to the canceler. A second 
sensor receives a noise n x which is correlated in some way with the noise w 0 . This sensor 
output provides the "reference input" to the canceler. The noise n x is filtered to produce 
an output y that is a close replica of n 0 . This output is subtracted from the primary 
input s + n 0 to produce the system output s + /?„ — y. 
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Figure 7. Adaptive Noise Canceler 



In the s)stem shown in Figure 7, the reference input is processed by our two 
dimensional least mean square filter. Thus the filter operates under changing conditions 
and will readjust itself continuously to minimize the error signal. 

The computer program which simulates this noise canceling model is provided 
in Appendix B. We will discuss the results and conclusions concerning various simu- 
lations in Chapter 4. 

3. Adaptive Line Enhancer 

A special case of adaptive noise canceling is when there is only one signal x(n) 
available which is contaminated by noise. In such a case, the signal .*(/;) provides its 
own reference signal >•(/?), which is taken to be a delayed replica of 
*(«): >•(/;) = x(n — A), as shown in Figure 8 on page 19. The adaptive filter will respond 
by canceling any components of the main signal that are in any way correlated with 
the secondary signal y(n) = x(n — A). Suppose the signal .*(/?) consists of two compo- 
nents: a narrowband component that has long-range correlations and a broadband 
component which will tend to have short-range correlations. One of these could repre- 
sent the desired signal and the other an undesired interfering noise. Suppose that the 
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Figure 8. Adaptive Line Enhancer 



delay A is selected so that it falls between the correlation lengths. Since A is longer than 
the effective correlation length of the broadband component, the delayed replica will be 
entirely uncorrclated with the broadband part of the main signal. The adaptive filter 
will not be able to respond to this component. On the other hand, since A is shorter 
than the correlation length of the narrowband component, the delayed replica that ap- 
pears in the secondary input will be correlated with narrowband part of the main signal 
and the filter will respond to cancel it. Note that if A is selected to be longer than both 
correlation lengths, the secondary input will become uncorrelated with the primary input 
and the adaptive filter will turn itself off. In the opposite case, when the delay A is se- 
lected to be less than both correlation lengths, then both components of the secondary 
signal will be correlated with the primary signal, and therefore the adaptive filter will 
respond to cancel the primary' x(n) completely. The computational algorithm for the 
adaptive line enhancer is shown in the following three equations: 



M M 

•*('0 = - "0 = - n\ - A) 

//i=U f?i— 0 



( 2 . 40 ) 
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e(n) = .v(/7) — .£(«) (2.41) 

h m (n + 1) = h m (n) + 2 n e(n) x(n — m — A) m = 0,1,2 ,. ,.,M (2-42) 

For this model we also developed a computer simulation which incorporates our 
two dimensional LMS algorithm and it is provided in Appendix C. The results and 
conclusions will again be discussed in Chapter 4. 
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III. TWO DIMENSIONAL ADAPTIVE RECURSIVE LEAST SQUARES 

ALGORITHM 



In the previous chapter, under the LMS algorithm, the available data samples were 
used in order to attempt to move the current estimate of the impulse response to the 
optimum value. This approach has the advantage of being simple to implement but 
carries with it the disadvantages that it can be slow to approach the optimal weight 
vector and, once close to it, will usually fluctuate around the optimal vector rather than 
actually converge to it due to the effects of approximations made in the estimate of the 
performance function gradient. 

To overcome these difficulties, we examine another approach in this chapter which 
uses the input data in such a way as to ensure optimality at each step. This alternative 
algorithm is based on the exact minimization of the least square criteria. This algorithm 
is known as recursive least squares (RLS). 

A. ONE DIMENSIONAL RECURSIVE LEAST SQUARES 

As in the LMS case, the one dimensional RLS algorithms is developed in several 
different ways [Refs. 5, 8 ]. Orfanidis [Ref. 6] provides a derivation which we will con- 
sider prior to our extension to two dimensions. The tapped delay line shown in 
Figure 9 on page 22 will provide the reference for the following discussion. We begin 
by replacing the LMS estimation criteria of MSE = £[<? 2 ] by 



n 




k = 0, !,...,« 



(3.2) 



where 



e{k) = x(k) — x(k) 



(3.3) 



and x(k) is the estimate of x(k) produced by the Mth-order Wiener Filter 




(3.4) 
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x(k) 



Figure 9. One Dimensional RLS Reference (TDL) 

Substituting equation (3.2) into (3.1) and setting the derivative with respect to h to zero 
we Find the least-squares analogs of the orthogonality equations 
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dMSE 

dh 



-2 




e(*)y(*) = o 



which may be rewritten in their normal equation form as follows 



(3.5) 



n 

i 

k = 0 



(x(k) - h T y(k))y(k) = 0 



E -. n 

y(k)y(k) T h = y W ) 
-I fe 0 



(3.6) 



(3.7) 



Defining the quantities 



R(n) = 


T]y(*) y( k ) T 


(3.8) 








r (") = 


n 

: y (*) 


(3.9) 



we can then write the normal equation as R(n) h = r(«), with solution h = R(n)- 1 r(«). 
Note that the n-dependence of R(n) and r(w) depend on the current time n, therefore, 



h(«) = R(rt) _1 r(«) = P(n)r(n) 



(3.10) 



where P(n ) = R(n)~'. These are the least squares analogs of the ordinary Wiener sol- 
ution, with R(n) and r(«) playing the role of the covariance matrix £[y(«)y(«) T ] and 
cross-correlation vector £[jr(«)y(n)], respectively. The RLS algorithm is obtained by 
writing equation (3.10) recursively in n and then using the following matrix inversion 
lemma [Ref. 5] 



(. A + BCD ) -1 = A~ ] - A~'B(DA~ l B + C~ h )~ ] DA~ ] 



we get the update equation for the P matrix 



P(n) = P(n- 1) - 



P(n - \)y(n)y(n) r P(n - 1) 
1 + y{n) T P{n- l)y(w) 



(3.11) 



(3.12) 
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Using the quantities in equations (3.8) and (3.9) to satisfy the recursions yields 



R(n) = R(n- 1) + y(n)y(«) r 



(3.13) 



r («) = r(* - 1) + *(rt)y(rt) 



(3.14) 



and substituting equations (3.12) and (3.14) into (3.10) after some mathematical manip- 
ulations we find 



which differs from the LMS alogorithm by the presence of the P(n), vice /i, in front of 
the weight correction term. Since P{n) = P{n)~ { is a measure of the covariance matrix 
£[y(«)y(")n, the presence of R(n)~ l makes the RLS algorithm behave like Newton's 
method, and hence has fast convergence properties. 

B. TWO DIMENSIONAL RECURSIVE LEAST SQUARES 

Although it is discussed briefly by Wellstead and Caldas Pintos [Ref. 14], limited 
information is available in the open literature in this area. In order to develop a two 
dimensional recursive least square (RLS) algorithm we will extend the one dimensional 
adaptive algorithm developed in the previous section in a method similar to that used for 
the LMS algorithm. 

As in Chapter 2, using Figure 10 on page 25 as a reference we see that the basic 
filter has two input images. The primary two dimensional input array, D, is the ideal 
image plus additive noise. The reference image X is the noise array. Each array is of 
dimension M by M. The filter mask, W, is N by N. 

The one dimensional RLS algorithm, equation (3.15), is the same as the one di- 
mensional LMS algorithm with the exception that P(n) replaces n. Therefore, if a 
method of developing a P matrix for the two dimensional case can be devised we can 
then use an algorithm similar to the two dimensional LMS algorithm to update the filter 
coefficients. First, let 



h (n) = h (n - 1) + P{n) y («) e{n) 



(3.15) 
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X ARRAY 




Figure 10. Two Dimensional RLS Reference 
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PJ. 0,0) p/0,1) 




Pj(0,n 2 -\) 



(3.16) 



Pj(N 2 - 1,0) P/N 2 -1,1) P/N 2 -1,A' 2 —1 ) 



which represents the P matrix on the jth iteration where j is defined as j = mM + n. 
As in the two dimensional LMS algorithm, we let the filter coefficient matrix be given 
as 



ff}(0,0) ^(0,1) 

wjtm WJL 1 , 1 ) 



Wj{0,N- 1) 
WjiUN-l) 



(3-17) 



lVj(N - 1 ,0) \Vj{N -1,1) ... Wj(N — \,N— 1) 



and the input data window as 



Xj = 



X(m,n ) X(m,n — 1) 



X{m,n- N+ 1) 



(3.18) 



X(m-N+l,n) X(m-N+\,n-l) ... X(m — N + \,n — N+ 1) 



In order to establish an update equation for the P matrix in which each element will 
have the proper dimensions, we use equations (3.17) and (3.18) and make the following 
transformations; the filter mask, W p equation (3.17) is transformed into a vector defined 
as 
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^( 0 , 0 ) 

^•(0,1) 



Wj(0,N- 1) 

^■( 1 , 0 ) 



WWj = 



^•(1,/V-l) 



(3.19) 



W^N — 1 , 0 ) 

»0(iV-U) 



Wj(N - \ ,N— 1) 



and we transform the Xj matrix, equation (3.18) into 



X{m,n) 
X{m,n — 1) 



X(m,n-N+ 1) 



XX j = 



X(m - N+ 1 ,«) 
X(m — N + 1 ,n — 1) 



X(m-N+ \,n-N+ 1) 



(3.20) 



Using the quanitities defined above, we can then write the P matrix update equation as 
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P h i(XXj)(XXj)P j _ i 
1 + (XXj)P h ,(XXj) 



(3.21) 



Pj = Pj-, ~ 



As discussed earlier, the one dimensional RLS weight updating algorithm difTers 
from the LMS in that it contains the P matrix while the LMS algorithm contains n. This 
must be considered in two dimensions since P is a matrix and /j. is a scalar. Therefore 
using the previous equations, we can obtain the update to the filter coefficients via 

WW m = WWj - (Pj)(ej)(XXj) (3.22) 



where e, is given as 



ej = D(m,n) - y{m,n ) 



(3.23) 



y{m,ri) is the convolution sum of the image in the reference input with the filter window 



y{m,n) 



A-l N-l 




WUJk) xn,k) 



(3.24) 



and the value P s is updated by equation (3.21). 

The cost function of the RLS criterion could be modified to include a windowing 
of the input data as follows 



USE 




e\k) 



This modification results in the following modified P matrix update: 




p^xxyxxjv^ 

1 + (XXj T )Pj_,(XXj) 



(3.25) 



(3.26) 



where q, the averaging or "forgetting factor" is a positive constant. It is usually chosen 
to be slightly less than one, thereby diminishing the contribution of the "older" data. 
The problem of data nonstationarity is the main reason for introducing this type of 
weighting factor.fRef. 6] It should be noted that the usual RLS algorithm is attained 
when q = 1 and that no changes in the amount of computation is required. 
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C. IMPLEMENTATION 
1. System Identification 

Following the mathematical development of the two dimensional recursive least 
squares algorithm, we implemented and tested it as was previously done with the two 
dimensional least mean square algorithm. Our initial testing was done in a system 
identification model as shown in Figure 11. As in the two dimensional LMS case, we 
used a known filter with the following output equation 

d{m,n) = .4 lV(m,n) + .61V(m-\,n) 

- .3 ir(m,n- 1) + .2 IV{m- \,n- 1) ' 




error 



Figure 11. Two Dimensional RLS System Identification 
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The two dimensional adaptive recursive least square filter output was given as 

y{m,n) = AO W{m,n ) + A 1 IV (m — 1 ,n) 

+ A2 IV(m t n — 1) + A3 W{m — 1,m — 1) < 3 ' 28 > 

The error signal was generated by taking the difference between the desired (known) 
signal, d(m,n), and the filter output, y{m,n). By using equations (3.21) and (3.22) to up- 
date the P matrix and the filter coefficients, we caused the filter weights to move toward 
their optimum values and the error signal to approach zero. 

The computer program for this model is provided is Appendix D. For this sys- 
tem, a 32x32 white gaussian noise matrix was generated for the input to both the known 
filter and the RLS adaptive filter. The rate of convergence is shown in Figure 12 on 
page 31 and Figure 13 on page 32 and as can be seen after approximately 70 iterations 
each of the coefficients had converged to with 10 -3 of the actual value and the error was 
10 - 5 . 

2. Adaptive Noise Canceler/ Adaptive Line Enhancer 

These two systems were both discussed in Chapter 2 with regard to implemen- 
tation of the two dimensional LMS algorithm. In order to further test the RLS algo- 
rithm we also implemented it within the noise canceler and the adaptive line enhancer. 
The computer programs wiiich w : ere used for the simulation are given Appendix E and 
Appendix F, respectively. The simulation results and discussion will be provided in 
Chapter 4. 
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Figure 12. RLS System Identification Rate of Convergence 
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Figure 13. RLS System Identification Rate of Convergence 
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IV. RESULTS AND CONCLUSIONS 



In this chapter, we present the experimental results from our implementation of the 
two dimensional LMS algorithm derived in Chapter II and the two dimensional RLS 
algorithm derived in Chapter HI. The two algorithms were used in a noise caneeler and 
an adaptive line enhancer, as discussed in Chapter II, to solve an image restoration 



problem. 

Figure 14 shows the original image of a house which is comprised of 128 by 128 
eight bit pixels. The image has a mean value of 131.68 and a variance of 3194.4. In 
Figure 15 on page 34, we have corrupted the original image with additive white gaussian 
noise (zero mean and variance 1600). 




Figure 14. Original Image 




Figure 15. Image plus Additive Noise 



A. NOISE CANCELER RESULTS 

In order to solve our image restoration problem, we initially implemented both the 
LMS and the RLS algorithms in. a noise canceler using a two by two filter matrix and 
a value of 70x1 0~ 9 for jx in the LMS algorithm, 

Figure 16 on page 35 shows the output when the LMS algorithm is processed 
through the image for one pass: this results in 16,384 updates of the filter corresponding 
to the 16,384 pixels. Figure 17 on page 36 is the results following two passe? through 
the image. No significant improvement was realized with more than two passes. 

Implementing a two by two RLS adaptive filter in the noise canceler produced very 
favorable results after one pass through the image. These results are shown in 
Figure IS on page 36 and compare well with the original image and the LMS output 
after two passes. 

For the LMS algorithm our best results were achieved with a value of 70x1 0~ 9 for 
(x, however in order to demonstrate the effect of changing this value Figure 19 on page 
37 represents the output with ,u ~ 35x10'“* and Figure 20 on page 37 is produced by a 
spatially varying normalized fx as discussed in Chapter II. It can be seen that for various 
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values of a, diiTerent features within the image were restored at different levels, in Fig- 
ure 16 on page 35, the edges and more detailed segments of the image were restored 
better than in Figure 19 on page 37, however areas of similar contrast were not as well 
restored. When using the normalized p the mean value more closely approached the 
original mean, however the variance was reduced. 

Finally, increasing the number of coefficients, using a three by three filter matrix in 
the LMS algorithm we obtained the results shown in Figure 21 on page 38 in one pass. 
This filter showed no significant improvement in the variance compared to the two by 
two filter; however, it did improve with respect to the mean. 

Table 1 on page 38 shows the restoration results comparing the mean and variance 
of the various outputs with that of the original image and the image plus noise. 




Figure 16. Restored Image: Noise Canceier/LMS (One pass) 
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Figure 37. Restored Image: Noise Canceler/LMS (Two passes) 




Figure 18. Restored Image: Noise Canceler/RLS (One pass) 




Figure 19. Restored image: Noise Canceier/LMS (different mu) 




Figure 20. Restored Image: Noise Canceler/LMS (normalized mu) 
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Figure 21. Restored Image: Noise Canceler/LMS (3 by 3) 



Table 1. NOISE CANCELER RESULTS 



TYPE OF FILTER 


MEAN 


VARIANCE 


Original Image 


131.68 


3194.4 


Image plus Noise 


134.25 


1647.7 


One pass 2x2 LMS 


114.52 


2034.3 


Two pass 2x2 LMS 


136.78 


3711.8 


One pass 2x2 RLS 


139.04 


3607.7 


One pass LMS(different 
mu) 


120.12 


1940.4 


One pass 

LMS(normalized mu) 


126.53 


1500.7 


One pass 3x3 LMS 


119.13 


2017.3 
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B, ADAPTIVE LINE ENHANCER RESULTS 



We next implemented our two algorithms in an adaptive line enhancer. In order to 
obtain a reference input from the primary input, we used a two dimensional delay oper- 
ator of (m-l,n). As in the noise canceler we again used a two by two filter matrix and 
this produced the results shown in Figure 22 on page 39 and Figure 23 on page 40 for 
the one pass LMS and the one pass RLS, respectively. No significant improvement was 
noted in the two pass LMS compared to the one pass LMS. 

Table 2 on page 40 shows the comparison of each output mean and variance with 
the mean and variance of the original image and the image plus noise. 




Figure 22. Restored Image: Adaptive Line Enhancer/LMS (One pass) 



1 




Table 2. ADAPTIVE LINE ENHANCER RESULTS 



TYPE OF FILTER j 


| MEAN 


VARIANCE 


Original Image 


! 131.68 1 


3 i 94.4 


Image plus Noise 


13425 


1(347.7 


One pass 2x2 LMS j 


f 125.03 


2521.8 


Two pass 2x2 LMS | 


1 126.94 j 

. _ . - - J 


2529.9 


One pass 2x2 RLS 


139.69 


2124.2 



C CONCLUSIONS 

In this thesis we have extended a one dimensional LMS and a one dimensional RLS 
adaptive algorithm to two dimensions. As we examine the results we see that many of 
the one dimensional comparisons between LMS and RLS are applicable to two dimen- 
sions: 

1. the rate of convergence of the RLS algorithm is in general faster than that of the 
LMS algorithm by an order of magnitude 

2, this superior performance of the RLS algorithm is attained at the expense of a large 
increase in the computational complexity 
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3. there are no approximations within the derivation of the RLS algorithm, therefore 
as the number of iterations approach infinity, the least squares estimate approaches 
the optimum Wiener value 

4. in the RLS algorithm the correction which is applied is computed using all past 
data whereas the LMS algorithm uses only the instantaneous sample and the error 
signal; this is not necessarily an advantage in image processing. 

The objectives of this thesis were successfully completed. Some suggestions for fu- 
ture work include: (i) to examine these two dimensional algorithms in other applica- 
tions, applying them to areas other than image processing, (ii) to extend the two 
dimensional LMS and RLS algorithms to n-dimensions and implement them, (iii) to 
analyze the extension of these one dimensional algorithms to two dimensions in the fre- 
quency domain. Multidimensional digital signal processing is being applied in many 
areas today, however the potential for future growth and applications appears unlimited. 
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APPENDIX A. LMS SYSTEM IDENTIFICATION PROGRAM 



SYSTEM IDENTIFICATION OF AN ADAPTIVE TWO DIMENSIONAL 
LMS FILTER VS. A KNOWN TWO DIMENSIONAL FILTER 

REAL MU 

DIMENSION A0( 0: 32,0: 32),A1(0: 32,0: 32),A2(0: 32,0: 32), 

* A3( 0: 32,0: 32),E(0: 32,0: 32),D(0: 32,0: 32),W(-1: 32,-1: 32), 

* Y(0: 32,0: 32) 

****** v^R I AB LE DEFINITIONS****** 

A=NOISE AMPLITUDE 
S=NOISE VARIANCE 
AM=NOISE MEAN 
IX=SEED 

A0 , A1 , A2 , A3=FILTER COEFFICIENTS 

MU=SCALING FACTOR WITHIN THE LMS ALGORITHM 

W=WHITE GAUSSIAN NOISE GENERATED BY SUBROUTINE 

Y=FILTER OUTPUT 

D=DESIRED OUTPUT 

E=ERROR(D - Y) 



******* INITIAL VA LUE S******* 

A=l. 

MU=. 01 
A0(0,0)=0. 

Al( 0, 0)=0. 

A2( 0,0)=0. 

A3(0,0)=0. 

S=l. 

AM=0. 

IX=65539 

*****OPEN FILE AND FILL INPUT MATRIX WITH WHITE 
GAUSSIAN NOISE******** 

OPEN (UNIT=80,FILE=' SYSID2 DATA ', STATUS= ’ NEW ’ ) 

DO 10 M=-l , 31 
DO 20 N=- 1,31 

IF( (M. LT. 0) . OR. (N. LT. 0) )THEN 
W(M,N)=0. 

ELSE 

CALL GAUSS( IX, S , AM,V) 

W(M,N)=V*A 

ENDIF 

CONTINUE 

CONTINUE 

******COMPUTE THE UPDATED FILTER COEFFICIENTS 

USING THE TWO DIMENSIONAL LMS ALGORITHM***** 

DO 40 K=0 , 31 
DO 50 J=0 ,31 

D(K, J)=. 4*W(K, J)+. 6*W(K-1 , J) 3*W(K,J-1)+. 2*W(K-1,J-1) 
Y(K, J)=A0(K, J)*W(K, J)+A1(K, J)*W(K-1 , J)+A2(K, J)*W(K, J-l)+ 
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* A3(K,J)*W(K-1,J-1) 

E(K, J)=D(K, J) -Y(K, J) 

A0(K,J+1)=A0(K,J)+MU*E(K,J)*W(K,J) 

A1(K,J+1)=A1(K,J)+MU*E(K,J)*W(K-1,J) 

A2(K, J+1)=A2(K, J)+MU*E(K, J)*W(K, J-l) 

A3(K, J+1)=A3(K, J)+MU*E(K, J)*W(K-1 , J-l) 

50 CONTINUE 

A0(K+1,0)=A0(K,J) 

A1(K+1,0)=A1(K,J) 

A2(K+1,0)=A2(K,J) 

A3(K+1,0)=A3(K, J) 

40 CONTINUE 

******0UTPUT RESULTS******* 

L=0 

DO 60 K=0 , 3 1 
DO 70 J=0,31 
L=L+1 

PRINT 15, L,E(K,J),A0(K,J),A1(K,J),A2(K,J),A3(K,J) 
WRITE (80,15) L,E(K,J),A0(K,J),A1(K,J),A2(K,J),A3(K,J) 
15 FORMAT (' ' , 1X,I4,2X,F8. 5,2X,F8. 5,2X,F8. 5,2X,F8. 5, 

* 2X,F8. 5) 

70 CONTINUE 

60 CONTINUE 
STOP 
END 
C 

SUBROUTINE GAUSS( IX, S , AM, V) 

A=0. 0 

DO 50 1=1,12 
CALL RANDU( IX, IY, Y) 

IX=IY 
50 A=A+Y 

V=(A-6. 0)*S+AM 
RETURN 
END 
C 
C 

SUBROUTINE RANDU( IX, IY,YFL) 

IY=IX*65539 

IF(IY)5,6,6 

5 I Y=IY+2 147483647+1 

6 YFL=IY 

YFL=YFL*. 4656613E-9 

RETURN 

END 
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APPENDIX B. LMS ALGORITHM IMPLEMENTED IN A NOISE 

CANCELER 



THIS IS A VAX/VMS FORTRAN PROGRAM THAT IMPLEMENTS A TWO 
DIMENSIONAL 2X2 ADAPTIVE LMS FILTER WITHIN A NOISE 
CANCELER 

INTEGER M,N,K,J 

BYTE B( 0: 127) ,BYTEE(0: 127) ,BYTEU(0: 127), 

INTEGERS INTE( 0: 127,0: 127),INTU(0: 127,0: 127), 

* IE , IU 

REAL*4 MU, A,FMINE ,FMAXE ,FMINU,FMAXU 
REAL*4 AA( 0: 3),E(0: 127,0: 127), U(0: 127,0: 127) 

REAL*4 W(-l: 127,-1: 127),Y(0: 127,0: 127),IM(0: 127,0: 127) 

****** VAR I AB LE DEFINITIONS****** 

A=N0ISE AMPLITUDE 
S=N0ISE VARIANCE 
AM=N0ISE MEAN 
AA=FILTER COEFFICIENTS 
IX=SEED 

MU=SCALING FACTOR WITHIN THE LMS ALGORITHM 
IM=INPUT IMAGE 

W=WHITE GAUSSIAN NOISE GENERATED BY SUBROUTINE 
U=SIGNAL PLUS N0ISE(IM + W) 

Y=FILTER OUTPUT 
E=ERROR( U - Y) 

FMINE , FMAXE , FMINU , FAMXU=PARAMETERS 

TO BE USED TO CONVERT DECIMAL DATA TO BYTE DATA 



******* initial va lue s ********* 

A=l. 

MU=35. E-9 
AA(0)=0. 

AA( 1)=0. 

AA(2)=0. 

AA( 3)=0. 

S=40. 

AM=0. 

IX=65539 
FMINE=1. E+10 
FMAXE=-1. E+10 
FMINU=1. E+10 
FMAXU=- 1. E+10 



******OPEN AN IMAGE FILE, CONVERT THE BYTE DATA INTO INTEGERS 
AND THEN PLACE THESE VALUES IN A MATRIX******* 

OPEN (UNIT=1 , NAME =' H0US1G. DAT* , TYPE ='0LD', ACCESS = 
’DIRECT', RECORDS IZE=32 , MAXREC=128) 

DO 100 M=0 , 127 

READ( 1 ' M+l) (B( J) , J=0, 127) 

DO 110 N=0 . 127 
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IF(B(N). LT. 0)THEN 
IM(M,N)=B(N)+256 

FT QF 

IM(M,N)=B(N) 

ENDIF 

110 CONTINUE 

100 CONTINUE 
C 

c ******* ADD WHITE GAUSSIAN NOISE TO THE IMAGE AND SET THE 

C VALUES OUTSIDE THE IMAGE TO ZERO******* 

DO 10 M=-l , 127 
DO 20 N=- 1,127 

IF ((M. LT. 0).OR. (N. LT. 0))THEN 
W(M,N)=0. 

IM(M,N)=0. 

ELSE 

CALL GAUSS ( IX, S , AM, V) 

W(M,N)=V*A 

ENDIF 

U(M,N)=IM(M,N)+W(M,N) 

20 CONTINUE 

10 CONTINUE 

C 

C ******USE THE LMS ADAPTIVE ALGORITHM TO UPDATE 

C THE FILTER COEFICIENTS********** 

DO 40 K=0 ,127 
DO 50 J=0 , 127 

Y(K, J)=AA( 0)*W(K, J)+AA( 1)*W(K-1 , J)+ 

* AA( 2)*W( K, J- 1)+AA( 3)*W(K-1 , J-l) 
E(K,J)=U(K,J)-Y(K,J) 

AA( 0 )=AA( 0 ) +MU*E(K , J)*W( K , J) 

AA(1)=AA( 1)+MU*E(K, J)*W(K-1, J) 

AA( 2)=AA( 2)+MU*E(K, J)*W(K, J-l) 

AA(3)=AA(3)+MU*E(K, J)*W(K-1, J-l) 

50 CONTINUE 

40 CONTINUE 

C 

C ********CHANGE SIGNAL PLUS NOISE AND ERROR 

C OUTPUT INTO BYTE DATA AND THEN WRITE 

TO A FILE* ******* 

OPEN (UNIT=2,NAME=' ERROR. DAT' ,TYPE='NEW' ,ACCESS= 

* 'DIRECT' , RECORDS I ZE=32,MAXREC=128) 

OPEN (UNIT=3,NAME=' SIGNOISE.DAT' ,TYPE='NEW' ,ACCESS= 

* 'DIRECT' ,RECORDSIZE=32,MAXREC=128) 

DO 500 1=0,127 

DO 550 J=0 , 127 

IF(E( I , J). LT. FMINE)THEN 
FMINE=E( I , J) 

ENDIF 

IF(E( I , J) . GT. FMAXE)THEN 
FMAXE=E( I , J) 

ENDIF 

550 CONTINUE 

500 CONTINUE 

IF (FMINE. LT. 0) THEN 
FMAXE=FMAXE - FM I NE 
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650 

600 



750 

700 



450 

400 

C 



ENDIF 

DO 600 1=0,127 
DO 650 J=0, 127 

IF(FMINE. LT. 0)THEN 
E( I , J)=E( I , J) -FMINE 
E( I , J)=E( I , J)*255. /FMAXE 
ELSE 

E( I , J)=(E( I , J)=FMINE)*255. /FMAXE 
ENDIF 

IE=NINT(E( I , J) ) 

IF(IE. GT. 127) THEN 
IE=IE-256 
ENDIF 

BYTEE( J) -IE 
CONTINUE 

WRITE ( 2 ' 1+1) (BYTEE(N) ,N=0 ,127) 
CONTINUE 
DO 700 1=0,127 
DO 750 J=0 , 127 

IF(U( I , J). LT. FMINU)THEN 
FMINU=U( I , J) 

ENDIF 

IF(U( I , J). GT. FMAXU)THEN 
FMAXU=U( I , J) 

ENDIF 

CONTINUE 

CONTINUE 

IF (FMINU. LT. 0) THEN 
FM AXU=FMAXU - FM I NU 
ENDIF 

DO 400 1=0,127 
DO 450 J=0 , 127 

IF(FMINU. LT. 0)THEN 
U( I , J)=U( I , J) -FMINU 
U(I,J)=U(I,J)*255. /FMAXU 
ELSE 

U( I , J)=(U( I , J)=FMINU)*255. /FMAXU 
ENDIF 

IU=NINT(U( I , J) ) 

IF( IU. GT. 127) THEN 
IU=IU-256 
ENDIF 

BYTEU( J) -IU 
CONTINUE 

WRITE ( 3 ' 1+1) (BYTEU(N) ,N=0, 127) 
CONTINUE 
CLOSE (UNIT=2) 

CLOSE (UNIT=3) 

STOP 

END 

SUBROUTINE GAUSS( IX,S,AM,V) 

AMP=0 . 0 
DO 50 1=1,12 

RANDOM=RAN( IX) 
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AMP=AMP+RANDOM 

CONTINUE 

V=( AMP-6. 0)*S+AM 

RETURN 

END 



ooo oooooooooooooooooo oooo 



APPENDIX C. LMS ALGORITHM IMPLEMENTED IN AN ADAPTIVE 

LINE ENHANCER 

THIS IS A VAX/ VMS FORTRAN PROGRAM THAT IMPLEMENTS A TWO 
DIMENSIONAL 2X2 ADAPTIVE LMS FILTER WITHIN AN ADAPTIVE 
LINE ENHANCER 

INTEGER M,N,K,J 

BYTE B( 0: 127) ,BYTEY( 0: 127) ,BYTEU(0: 127) 

INTEGER* 4 INTY(0: 127,0: 127),INTU(0: 127,0: 127),IY,IU 
REAL*4 MU , A , FMINY , FMAXY , FMINU , FMAXU 
REAL*4 AA( 0: 3),E(0: 127,0: 127), U(0: 127,0: 127) 

REAL*4 W( - 1: 127,-1: 127),Y(0: 127,0: 127),IM(0: 127,0: 127) 

******VARI AB LE DEFINITI ONS****** 

A=NOISE AMPLITUDE 
S=NOISE VARIANCE 
AM=NOISE MEAN 
AA=FILTER COEFFICIENTS 
IX=SEED 

MU=SCALING FACTOR WITHIN THE LMS ALGORITHM 
I M= INPUT IMAGE 

W=WHITE GAUSSIAN NOISE GENERATED BY SUBROUTINE 
U=SIGNAL PLUS NOISE(IM + W) 

Y=FILTER OUTPUT 

WW=DELAYED VERSION OF SIGNAL PLUS NOISE(U) 

E=ERROR(U - Y) 

FMINY , FMAXY , FMINU , F AMXU=P ARAMETERS 

TO BE USED TO CONVERT DECIMAL DATA TO BYTE DATA 



******* INI T I AL VALUES********* 
A=l. 

MU=35. E-9 
AA(0)=0. 

AA( 1)=0. 

AA(2)=0. 

AA(3)=0. 

S=40. 

AM=0. 

IX=65539 
FMINU=1. E+10 
FMAXU=-1. E+10 
FMINY=1. E+10 
FMAXY=-1. E+10 



******0PEN AN IMAGE FILE, CONVERT THE BYTE DATA INTO INTEGERS 
AND THEN PLACE THESE VALUES IN A MATRIX******* 

OPEN (UNIT=1 , NAME ='H0US1G. DAT' , TYPE ='OLD', ACCESS = 

* 'DIRECT', RECORDSIZE=32, MAXREC=128) 

DO 100 M=0 , 127 

READ(1'M+1) (B( J) , J=0, 127) 

DO 110 N=0, 127 
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IF(B(N). LT. 0)THEN 
IM(M,N)=B(N)+256 

VT QT7 

IM(M,N)=B(N) 

ENDIF 

110 CONTINUE 

100 CONTINUE 
C 

c *******ADD WHITE GAUSSIAN NOISE TO THE IMAGE AND SET THE 
C VALUES OUTSIDE THE IMAGE TO ZERO******* 

DO 10 M=-l , 127 
DO 20 N=- 1,127 

IF ( (M. LT. 0).OR. (N. LT. 0))THEN 
W(M,N)=0. 

IM(M,N)=0. 

ELSE 

CALL GAUSS( IX,S, AM,V) 

W(M,N)=V*A 

ENDIF 

U(M,N)=IM(M,N)+W(M,N) 

WW(M,N)=U(M-1 jN) 

20 CONTINUE 

10 CONTINUE 
C 

C ******USE THE LMS ADAPTIVE ALGORITHM TO UPDATE 

C THE FILTER COEFICIENTS********** 

DO 40 K=0, 127 ■ 

DO 50 J=0 , 127 

Y(K,J)=AA(0)*WW(K,J)+AA(1)*WW(K-1,J)+ 

* AA( 2) *WW(K , J- 1 )+AA( 3 ) *WW( K- 1 , J- 1 ) 
E(K,J)=U(K,J)-Y(KjJ) 

AA( 0 )=AA( 0 )+MU*E ( K , J)*WW( K , J ) 

AA( 1)=AA( 1)+MU*E(K, J)*WW(K-1 , J) 

AA( 2)=AA( 2)+MU*E(K, J)*WW(K, J-l) 

AA( 3)=AA( 3)+MU*E(K, J)*WW(K-1, J-l) 

50 CONTINUE 

40 CONTINUE 

C 

C ********CHANGE SIGNAL PLUS NOISE AND FILTER OUTPUT 

C INTO BYTE DATA AND THEN WRITE TO A FILE****** 

OPEN (UNIT=3,NAME=' SIGNOISE.DAT' ,TYPE='NEW' ,ACCESS= 

* 'DIRECT' ,RECORDSIZE=32 ,MAXREC=128) 

OPEN (UNIT=4,NAME=' FILTERED. DAT' ,TYPE='NEW' ,ACCESS= 

* 'DIRECT' ,RECORDSIZE=32 ,MAXREC=128) 

DO 700 1=0,127 

DO 750 J=0 , 127 

IF( U( I , J ) . LT. FMINU ) THEN 
FMINU=U( I , J) 

ENDIF 

IF(U( I , J). GT. FMAXU)THEN 
FMAXU=U( I , J) 

ENDIF 

750 CONTINUE 

700 CONTINUE 

IF (FMINU. LT. 0) THEN 
FMAXU=FMAXU -FMINU 
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450 

400 



850 

800 



950 

900 

C 



ENDIF 

DO 400 1=0,127 
DO 450 J=0 , 127 

IF( FMINU.LT. 0)THEN 
U( I , J)=U( I , J) -FMINU 
U( I , J)=U( I , J)*255. /FMAXU 
ELSE 

U( I , J)=(U( I , J)=FMINU)*255. /FMAXU 
ENDIF 

IU=NINT(U( I , J) ) 

IF( IU. GT. 127) THEN 
IU=IU-256 
ENDIF 

BYTEU( J) -IU 
CONTINUE 

WRITE ( 3 ' 1+1) (BYTEU(N) ,N=0, 127) 
CONTINUE 
DO 800 1=0,127 
DO 850 J=0 ,127 

IF( Y( I , J) . LT. FMINY)THEN 
FMINY=Y( I , J) 

ENDIF 

IF(Y( I , J). GT. FMAXY)THEN 
FMAXY=Y( I , J) 

ENDIF 

CONTINUE 

CONTINUE 

IF (FMINY. LT. 0) THEN 
FMAXY=FMAXY - FM I NY 
ENDIF 

DO 900 1=0,127 
DO 950 J=0 , 127 

IF(FMINY. LT. 0)THEN 
Y(I,J)=Y(I,J) -FMINY 
Y( I , J)=Y( I , J)*255. /FMAXY 
ELSE 

Y( I , J)=(Y( I , J)=FMINY)*255. /FMAXY 
ENDIF 

IY=NINT( Y( I , J) ) 

IF(IY. GT. 127) THEN 
IY=IY-256 
ENDIF 

BYTEY( J) -IY 
CONTINUE 

WRITE (4'I+1) (BYTEY(N) ,N=0, 127) 
CONTINUE 
CLOSE (UNIT=3) 

CLOSE (UNIT=4) 

STOP 

END 

SUBROUTINE GAUSS ( IX, S , AM,V) 

AMP=0. 0 
DO 50 1=1,12 

RANDOM=RAN( IX) 
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AMP=AMP+RANDOM 

CONTINUE 

V=( AMP-6. 0)*S+AM 

RETURN 

END 



OO Ln OO OOOOOOOOOOOOOOOOO OOOO 



APPENDIX D. RLS SYSTEM IDENTIFICATION 

SYSTEM IDENTIFICATION OF AN ADAPTIVE RLS 2X2 FILTER VERSUS 
2X2 RLS FILTER VERSUS A A 2X2 FILTER WITH 
KNOWN COEFFICIENTS 



DIMENSION W( -1: 32,-1: 32),P(0: 3,0: 3),AA(0: 3),X(0: 3) 

****** VAR I AB LE DEFINITIONS****** 

A=N0ISE AMPLITUDE 
S=N0ISE VARIANCE 
AM=N0ISE MEAN 
IX=SEED 

AA=FILTER COEFFICIENTS 

W=WHITE GAUSSIAN NOISE GENERATED BY SUBROUTINE 
D=DESIRED OUTPUT 
Y=FILTER OUTPUT 
E=ERR0R( D - Y) 

MXM = ARRAY SIZE 
NXN = FILTER SIZE 
Q=WEIGHTING FACTOR 
P=INVERSE OF COVARIANCE MATRIX 



INPUT INITIAL VALUES ************ 

M=32 

N=2 

A=l. 

AA( 0)=0. 

AA( 1)=0. 

AA( 2)=0. 

AA(3)=0. 

Q=l. 

S=1 

AM=0 

IX=65539 

****** BUILD INITIAL 'P' MATRIX ********* 

DO 1 K=0,(N**2-1) 

DO 5 J=0,(N**2-1) 

IF(K. EQ. J)THEN 
P(K, J)=100. 

ELSE 

P(K, J)=0. 0 
END IF 
CONTINUE 
CONTINUE 

OPEN (UNIT=63 ,FILE=' RLSID2D2 DATA ', STATUS= ' OLD ' ) 

******* CALCULATE INPUT MATRIX (WHITE GAUSSIAN NOISE) ****** 
DO 10 K=-1,M-1 
DO 20 J=-1,M-1 

IF((K.LT. 0).OR. (J. LT. 0))THEN 
W(K, J)=0. 
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20 

10 

C 

C 

C 

c 

50 

40 

15 

C 

C 

C 

20 

30 

40 

50 

60 

70 

80 



ELSE 

' CALL GAUSS ( IX, S, AM, V) 

W(K,J)=V*A 

ENDIF 

CONTINUE 

CONTINUE 

********* CALCULATE ERROR BETWEEN ADAPTIVE FILTER OUTPUT 
AND KNOWN FILTER OUTPUT. USE THIS ERROR TO UPDATE 
FILTER COEFFICIENTS AND THE 'P' MATRIX ************* 
DO 40 K=0 , M-l 
DO 50 J=0,M-1 
L=K*M+J 

D= 4*W(K,J)+. 6*W(K-1,J)-. 3*W(K,J-1)+. 2*W(K-1,J-1) 
Y=AA(0)*W(K,J)+AA(1)*W(K-1,J)+AA(2)*W(K,J-1)+ 

* AA(3)*W(K-1,J-1) 

E=D-Y 

X(0)=W(K,J) 

X(1)=W(K-1,J) 

X(2)=W(K, J-l) 

X(3)=W(K-1,J-1) 

CALL RLS(P,X,AA,N,E,Q) 

PRINT 15, L,E,AA(0) ,AA( 1) , AA(2) ,AA(3) 

WRITE (63,15) L,E,AA(0) ,AA( 1) ,AA(2) ,AA(3) 

CONTINUE 

CONTINUE 

FORMAT (' ' ,2X,I3,2X,F9. 5,2X,F9. 5,2X,F9. 5,2X,F9. 5,2X,F9. 5) 

STOP 

END 



********* TO COMPUTE THE GAIN MATRIX ********** 
SUBROUTINE RLS(P,X,AA,N,E,Q) 

DIMENSION P(0: 3,0: 3),X(0: 3),AA(0: 3),TEMY(0: 3,0: 3), 
* GAMA(0: 3) ,TEMC(0: 3) 

DO 30 L=0 , 3 
TEM=0. 

DO 20 K=0 , 3 
TEM=TEM+X(K)*P(K,L) 

GAMA(L)=TEM 

GAM=Q 

DO 40 K=0 , 3 
GAM=GAM+GAMA( K ) *X( K ) 

DO 60 L=0 , 3 
TEM=0. 

DO 50 K=0 , 3 
TEM=TEM+P(L,K)*X(K) 

TEMC( L)=TEM 
DO 70 L=0 , 3 
DO 70 K=0 , 3 

TEMY(L,K)=TEMC(L)*GAMA(K) 

DO 80 K=0,3 
DO 80 L=0,3 

P(K,L)=(P(K,L)-(TEMY(K,L)/GAM))/Q 
DO 100 K=0,3 
TEM=0. 
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DO 90 L=0 , 3 



90 


TEM=TEM+P(K,L)*X(L) 


100 


TEMC(K)=TEM 




DO 110 K=0 , 3 


110 


AA(K)=AA(K)+TEMC(K)*E 




Q= 99*Q+. 01 




RETURN 




END 
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APPENDIX E. RLS ALGORITHM IMPLEMENTED IN A NOISE 

CANCELER 

THIS IS A VAX/ VMS FORTRAN PROGRAM THAT IMPLEMENTS A TWO 
DIMENSIONAL 2X2 ADAPTIVE RLS FILTER WITHIN A NOISE 
CANCELER 

INTEGER M,N,K, J 

BYTE B( 0: 127) ,BYTEE(0: 127) ,BYTEU(0: 127) 

INTEGERS INTEC 0: 127,0: 127),INTU(0: 127,0: 127),IE,IU 
REAL*4 A,FMINE ,FMAXE ,FMINU,FMAXU 
REAL*4 AA(0: 3) ,E(0: 127,0: 127), U( 0: 127,0: 127) 

REAL*4 P(0: 3,0: 3),X(0: 3) ,IM(0: 127,0: 127) 

REAL*4 WC-1: 127,-1: 127),Y(0: 127,0: 127) 

C 

c ******variable DEFINITIONS****** 

C A=NOISE AMPLITUDE 

C S=NOISE VARIANCE 

C AM=NOISE MEAN 

C AA=FILTER COEFFICIENTS 

C IX=SEED 

C I M= INPUT IMAGE 

C W=WHITE GAUSSIAN NOISE GENERATED BY SUBROUTINE 
C U=SIGNAL PLUS NOISE(IM + W) 

C Y=FILTER OUTPUT 

C E=ERROR( U - Y) 

C Q=WEIGHTING FACTOR 

C P=INVERSE OF COVARIANCE MATRIX 

C FMINE, FMAXE , FMINU , FAMXU=PARAMETERS 

C TO BE USED TO CONVERT DECIMAL DATA TO BYTE DATA 

C MXM=ARRAY SIZE 

C NXN=FILTER SIZE 

C 

c ******* initial values********* 

M=128 

N=2 

Q=l. 

A=l. 

AA( 0)=0. 

AA( 1)=0. 

AA( 2)=0. 

AA( 3)=0. 

S=40. 

AM=0. 

IX=65539 
FMINE=1. E+10 
FMAXE=-1. E+10 
FMINU=1. E+10 
FMAXU=- 1. E+10 
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c 

C ****** open AN IMAGE FILE, CONVERT THE BYTE DATA INTO INTEGERS 

C AND THEN PLACE THESE VALUES IN A MATRIX******* 

OPEN (UNIT=1, NAME =' HOUS1G. DAT* , TYPE ='OLD', ACCESS = 

* 'DIRECT', RECORDSIZE=32 , MAXREC=128) 

DO 100 K=0 , 127 

READ( 1 ' K+l) (B( L) , L=0 , 127) 

DO 110 J=0 , 127 

IF(B( J). LT. 0)THEN 
IM(K, J)=B( J)+256 
ELSE 

IM(K,J)=B(J) 

END IF 

110 CONTINUE 

100 CONTINUE 

CLOSE (UNIT=1) 

C 

C ****** BUILD INITIAL ' P ' MATRIX ********* 

DO 1 K=0 , (N**2- 1) 

DO 5 J=0 , ( N**2- 1) 

IF(K. EQ. J)THEN 
P(K, J)=100. 

ELSE 

P(K, J)=0. 0 
END IF 

5 CONTINUE 

1 CONTINUE 

C 

C ******* ADD WHITE GAUSSIAN NOISE TO THE IMAGE AND SET THE 

C VALUES OUTSIDE THE IMAGE TO ZERO******* 

DO 10 K= - 1 , M - 1 
DO 20 J=- 1 , M- 1 

IF ( (K. LT. 0). OR. (J. LT. 0))THEN 
W(K, J)=0. 

IM(K, J)=0. 

ELSE 

CALL GAUSS( IX, S , AM, V) 

W(K, J)=V*A 
ENDIF 

U(K, J)=IM(K, J)+W(K, J) 

20 CONTINUE 

10 CONTINUE 

C 

C ********* CALCULATE ERROR BETWEEN ADAPTIVE FILTER OUTPUT 

C AND KNOWN FILTER OUTPUT. USE THIS ERROR TO UPDATE 

C FILTER COEFFICIENTS AND THE ' P ' MATRIX ************* 

DO 40 K=0,M-1 
DO 50 J=0 , M- 1 
L=K*M+J 

Y=AA( 0)*W(K, J)+AA( 1)*W(K- 1 , J)+AA( 2)*W(K, J-l)+ 

* AA(3)*W(K-1, J-l) 

E(K, J)=U(K , J) -Y(K, J) 

EE=E(K, J) 

X(0)=W(K,J) 

X( 1)=W(K-1 , J) 

X(2)=W(K, J-l) 
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X(3)=W(K-1 , J-l) 

CALL RLS ( P , X , AA j N j EE j Q ) 

CONTINUE 

CONTINUE 

********change SIGNAL PLUS NOISE AND ERROR OUTPUT 
INTO BYTE DATA AND WRITE TO A FILE***** 

OPEN (UNIT=2,NAME=* ERROR. DAT’ ,TYPE='NEW' ,ACCESS= 

* 'DIRECT* ,REC0RDSIZE=32,MAXREC=128) 

OPEN (UNIT=3.NAME=' SIGNOISE. DAT* ,TYPE=’NEW* , ACCESS 

* 'DIRECT* ,RECORDSIZE=32,MAXREC=128) 

DO 500 1=0,127 

DO 550 J=0, 127 

IF(E( I , J). LT. FMINE)THEN 
FMINE=E( I , J) 

END IF 

IF(E(I ,J). GT. FMAXE)THEN 
FMAXE=E( I , J) 

END IF 
CONTINUE 
CONTINUE 

IF (FMINE. LT. 0) THEN 
FM AXE=FM AXE - F M I NE 
ENDIF 

DO 600 1=0,127 
DO 650 J=0 ,127 

IF(FMINE. LT. 0)THEN 
E( I , J)=E( I , J) -FMINE 
E( I , J)=E( I , J)*255. /FMAXE 
ELSE 

' E(I,J)=(E(I , J)=FMINE)*255. /FMAXE 
ENDIF 

IE=NINT(E( I , J) ) 

IF(IE. GT. 127) THEN 
IE=IE-256 
ENDIF 

BYTEE( J) -IE 
CONTINUE 

WRITE (2’I+1) (BYTEE(N) ,N=0, 127) 

CONTINUE 
DO 700 1=0,127 
DO 750 J=0, 127 

IF(U( I , J). LT. FMINU)THEN 
FMINU=U( I , J) 

ENDIF 

IF(U(I,J). GT. FMAXU)THEN 
FMAXU=U( I , J) 

ENDIF 

CONTINUE 

CONTINUE 

IF (FMINU. LT. 0) THEN 
FM AXU=FMAXU - FM I NU 
ENDIF 

DO 400 1=0,127 
DO 450 J=0 , 127 

IF( FMINU. LT. 0)THEN 
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U( I , J)=U( I , J) -FMINU 
U(I,J)=U(I,J)*255. /FMAXU 
ELSE 

* U( I , J)=(U( I , J)=FMINU)*255. /FMAXU 
ENDIF 

IU=NINT( U( I , J) ) 

IF(IU. GT. 127) THEN 
IU=IU-256 
ENDIF 

BYTEU( J) -IU 
CONTINUE 

WRITE ( 3 ' 1+1) (BYTEU(N) ,N=0 ,127) 

CONTINUE 
CLOSE (UNIT=2) 

CLOSE (UNIT=3) 

STOP 

END 

********* TO COMPUTE THE GAIN MATRIX ********** 
SUBROUTINE RLS( P, X , AA ,N ,E ,Q) 

DIMENSION P( 0: 3,0: 3) ,X(0: 3) ,AA(0: 3) ,TEMY(0: 3,0: 3) , 
* GAMA( 0:3) ,TEMC( 0:3) 

DO 30 L=0 , 3 
TEM=0. 

DO 20 K=0 , 3 
TEM=TEM+X(K)*P(K,L) 

GAMA(L)=TEM 

GAM=Q 

DO 40 K=0 , 3 

GAM=G AM+G AM A( K)*X(K) 

DO 60 L=0 , 3 
TEM=0. 

DO 50 K=0 , 3 
TEM=TEM+P(L,K)*X(K) 

TEMC(L)=TEM 
DO 70 L=0 , 3 
DO 70 K=0 , 3 

TEMY(L,K)=TEMC(L)*GAMA(K) 

DO 80 K=0 , 3 
DO 80 L=0 , 3 

P(K,L)=(P(K,L^-(TEMY(K,L)/GAM))/Q 
DO 100 K=0 , 3 
TEM=0. 

DO 90 L=0 , 3 
TEM=TEM+P(K ,L)*X(L) 

TEMC(K)=TEM 

DO 110 K=0 , 3 

AA(K)=AA(K)+TEMC(K)*E 

Q— . 99*Q+. 01 

RETURN 

END 
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APPENDIX F. RLS ALGORITHM IMPLEMENTED IN AN ADAPTIVE 



LINE ENHANCER 

THIS IS A VAX/VMS FORTRAN PROGRAM THAT IMPLEMENTS A TWO 


APP00010 


DIMENSIONAL 2X2 ADAPTIVE RLS FILTER WITHIN AN ADAPTIVE 


APP00020 


LINE ENHANCER 


APP00030 


INTEGER M,N,K, J 


APP00040 

APP00050 


BYTE B( 0: 127) ,BYTEU( 0: 127) ,BYTEY(0: 127) 


APP00060 


INTEGERS INTY( 0: 127,0: 127),INTU(0: 127,0: 127),IU,IY 


APP00070 


REAL*4 A,FMINU,FMAXU,FMINY,FMAXY 


APP00080 


REAL*4 AA( 0: 3) ,E(0: 127,0: 127) ,U(0: 127,0: 127) 


APP00090 


REAL*4 P( 0: 3,0: 3),X(0: 3),IM(0: 127,0: 127) 
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REAL* 4 W( -1: 127,-1: 127), Y(0: 127,0: 127) 
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******VARI ABLE DEFINITIONS****** 
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A=NOISE AMPLITUDE 


APP00140 


S=NOISE VARIANCE 
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AM=NOISE MEAN 


APP00160 


AA=FILTER COEFFICIENTS 


APP00170 


IX=SEED 


APP00180 


I M= INPUT IMAGE 
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W=WHITE GAUSSIAN NOISE GENERATED BY SUBROUTINE 
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U=SIGNAL PLUS NOISE(IM + W) 
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Y=FILTER OUTPUT 
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E=ERROR( U - Y) 
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WW=DELAYED VERSION OF SIGNAL PLUS NOISE(U) 
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Q=WEIGHTING FACTOR 


APP00250 


P=INVARSE OF COVARIANCE MATRIX 
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FMINU , FAMXU , FM INY , FMAXY=P ARAMETERS 
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TO BE USED TO CONVERT DECIMAL DATA TO BYTE DATA 
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MXM=ARRAY SIZE 
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NXN=FILTER SIZE 
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******* INITIAL VALUE S********* 
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M=128 
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N=2 


APP00340 


Q=l. 
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A=l. 
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AA( 0)=0. 


APP00370 


AA( 1)=0. 


APP00380 


AA(2)=0. 
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AA( 3)=0. 
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S=40. 


APP00410 


AM=0. 


APP00420 


IX=65539 


APP00430 


FMINU=1. E+10 


APP00440 


FMAXU=-1. E+10 


APP00450 


FMINY=1. E+10 


APP00460 


FMAXY=-1. E+10 
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C ****** OPEN AN IMAGE FILE, CONVERT THE BYTE DATA INTO INTEGERS 

C AND THEN PLACE THESE VALUES IN A MATRIX******* 

OPEN (UNIT=1, NAME =' HOUS1G. DAT* , TYPE ='OLD’ , ACCESS = 

* ’DIRECT’, RECORDS I ZE=32, MAXREC=128) 

DO 100 K=0 , 127 

READ( l’K+1) (B(L) ,L=0, 127) 

DO 110 J=0 , 127 

IF(B( J). LT. 0)THEN 
IM(K, J)=B( J)+256 
ELSE 

IM(K, J)=B( J) 

END IF 

110 CONTINUE 

100 CONTINUE 

CLOSE (UNIT=1) 

C 

C ****** BUILD INITIAL * P ' MATRIX ********* 

DO 1 K=0 , ( N**2 - 1 ) 

DO 5 J=0 , ( N**2 - 1 ) 

IF(K. EQ. J)THEN 
P(K,J)=100. 

ELSE 

P(K, J)=0. 0 
ENDIF 

5 CONTINUE 

1 CONTINUE 

C 

C ******* ADD WHITE GAUSSIAN NOISE TO THE IMAGE AND SET THE 

C VALUES OUTSIDE THE IMAGE TO ZERO******* 

DO 10 K=-2 , M- 1 
DO 20 J=-2,M-1 

IF ((K. LT. 0).OR. (J. LT. 0))THEN 
W(K, J)=0. 

IM(K,J)=0. 

ELSE 

CALL GAUSS ( IX, S , AM, V) 

W(K,J)=V*A 

ENDIF 

U(K, J)=IM(K, J)+W(K, J) 

20 CONTINUE 

10 CONTINUE 

C 

c *****coMPUTE THE DELAYED VERSION OF THE SIGNAL 

C PLUS NOISE MATR I X ( U ) ****** 

DO 45 K=- 1,127 
DO 46 J=-l, 127 

IF((M. LT. -l).OR. (N. LT. -1) THEN 
WW(K, J)=0. 

ELSE 

WW(K, J)=U(K-1 , J) 

ENDIF 

46 CONTINUE 

45 CONTINUE 

C 

C ********* CALCULATE ERROR BETWEEN ADAPTIVE FILTER OUTPUT 
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C AND KNOWN FILTER OUTPUT. USE THIS ERROR TO UPDATE 

C FILTER COEFFICIENTS AND THE * P* MATRIX ************* 

DO 40 K=0,M-1 
DO 50 J=0,M-1 

Y=AA(0)*WW(K,J)+AA(1)*WW(K-1,J)+AA(2)*WW(K,J-1)+ 

* AA(3)*WW(K-1,J-1) 

E(K,J)=U(K,J)-Y(K,J) 

EE=E(K, J) 

X(0)=WW(K,J) 

X(1)=WW(K-1,J) 

X(2)=WW(K, J-l) 

X(3)=WW(K-1, J-l) 

CALL RLS(P ,X, AA,N,EE ,Q) 

50 CONTINUE 

40 CONTINUE 

C 

C ********CHANGE SIGNAL PLUS NOISE AND FILTER OUTPUT 

C INTO BYTE DATA AND WRITE TO A FILE***** 

OPEN (UNIT=3 ? NAME=' SIGNOISE. DAT* ,TYPE='NEW' ,ACCESS= 

* 'DIRECT' ,RECORDSIZE=32 ,MAXREC=128) 

OPEN (UNIT=4.NAME=' FILTERED. DAT' ,TYPE='NEW' ,ACCESS= 

* 'DIRECT , RECORDSIZE=32 ,MAXREC=128) 

DO 700 1=0,127 

DO 750 J=0 , 127 

IF(U( I , J). LT. FMINU)THEN 
FMINU=U( I , J) 

ENDIF 

IF(U( I , J). GT. FMAXU)THEN 
FMAXU=U( I , J) 

ENDIF 

750 CONTINUE 

700 CONTINUE 

IF (FMINU. LT. 0) THEN 
FMAXU=FMAXU - FM I NU 
ENDIF 

DO 400 1=0,127 
DO 450 J=0 , 127 

IF( FMINU. LT. 0)THEN 
U( I , J)=U( I , J) -FMINU 
U( I , J)=U( I , J)*255. /FMAXU 
ELSE 

U( I , J)=(U( I , J)=FMINU)*255. /FMAXU 
ENDIF 

IU=NINT(U( I , J) ) 

IF( IU. GT. 127) THEN 
IU=IU-256 
ENDIF 

BYTEU( J) -IU 
450 CONTINUE 

WRITE ( 3 ' 1+1) ( BYTEU( N) , N=0 , 127 ) 

400 CONTINUE 

DO 800 1=0,127 
DO 850 J=0 , 127 

IF(Y( I , J) . LT. FMINY)THEN 
FMINY=Y( I , J) 

ENDIF 
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850 
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950 
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C 

C 

20 

30 

40 

50 

60 

70 

80 



IF(Y(I,J). GT. FMAXY)THEN 
FMAXY=Y( I , J) 

ENDIF 

CONTINUE 

CONTINUE 

IF (FMINY. LT. 0) THEN 
FMAXY=FMAXY-FMINY 
ENDIF 

DO 900 1=0,127 
DO 950 J=0 ,127 

IF(FMINY. LT. 0)THEN 
Y(I,J)=Y(I,J) -FMINY 
Y(I,J)=Y(I,J)*255. /FMAXY 
ELSE 

Y(I,J)=(Y(I,J)=FMINY)*255. /FMAXY 
ENDIF 

IY=NINT(Y( I , J) ) 

IF(IY.GT. 127) THEN 
IY=IY-256 
ENDIF 

BYTEY( J) -IY 
CONTINUE 

WRITE (4*1+1) (BYTEY(N) ,N=0, 127) 
CONTINUE 
CLOSE (UNIT=3) 

CLOSE (UNIT=4) 

STOP 

END 



********* to COMPUTE THE GAIN MATRIX ********** 
SUBROUTINE RLS(P,X, AA,N,E ,Q) 

DIMENSION P( 0: 3,0: 3) ,X(0: 3) ,AA(0: 3) ,TEMY(0: 3,0: 3) , 
GAMA(0: 3) ,TEMC(0: 3) 

DO 30 L=0 , 3 
TEM=0. 

DO 20 K=0 , 3 
TEM=TEM+X(K)*P(K,L) 

GAMA(L)=TEM 

GAM=Q 

DO 40 K=0 , 3 
GAM=GAM+GAMA( K)*X( K) 

DO 60 L=0 , 3 
TEM=0. 

DO 50 K=0 , 3 
TEM=TEM+P(L,K)*X(K) 

TEMC(L)=TEM 
DO 70 L=0,3 
DO 70 K=0, 3 

TEMY(L,K)=TEMC(L)*GAMA(K) 

DO 80 K=0,3 
DO 80 L=0 ,3 

P(K,L)=(P(K,L)-(TEMY(K,L)/GAM) )/Q 
DO 100 K=0 , 3 
TEM=0. 
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DO 90 L=0 , 3 
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90 


TEM=TEM+P(K,L)*X(L) 
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100 


TEMC(K)=TEM 
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DO 110 K=0,3 
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110 


AA(K)=AA(K)+TEMC(K)*E 


APP02180 




Q=. 99*Q+. 01 
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RETURN 
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END 


APP02210 



63 



LIST OF REFERENCES 



1. Andrews, H. and Hunt, B., Digital Image Restoration, Prentice-Hall, Inc., Reading, 
Massachusetts, 1987. 

2. Gonzalez, R. and Wintz, P., Digital Image Processing, Addison-Wesley Publishing 
Co., Englewood Cliffs, New Jersey, 1977. 

3. Rosenfeld, A. and Kak, A., Digital Picture Processing, Academic Press, New York, 
1976. 

4. Murphy, P. K., Survey of Image Restoration Techniques, The Johns Hopkins Uni- 
versity, Applied Physics Laboratory, Laurel, Maryland, July 1988. 

5. Treichler, J., Johnson, C., and Larimore, M., Theory and Design of Adpaptive Filters, 
John Wiley and Sons Inc., New York, 1984. 

6. Orfanidis, S. J., Optimum Signal Processing, an Introduction, Macmillian, 1985. 

7. Giordano, A. and Hsu, F., Least Square Estimation with Applications to Digital 
Signal Processing, John Wiley and Sons Inc., New York, 1985. 

8. Haykin, S., Introduction To Adaptive Filters, Macmillian Publishing Company, New 
York, 1984. 

9. Hadhoud, M. M. and Thomas D. W., "The Two-Dimensional Adaptive LMS Al- 
gorithm," IEEE Trans . Circuits Syst., Vol. CAS-35, no. 5, pp. 485-494, May 1988. 

10. Ekstrom, M. P., "Realizable Wiener Filtering in Two Dimensions," IEEE Trans. 
Acoust., Speech, Signal Processing, Vol. ASSP-30, no. 30, pp. 31-39, February 1982. 



64 



11. Florian, S. and Feuer, A., "Performance Analysis of the LMS Algorithm with a 
Tapped Delay Line (Two Dimensional Case)," IEEE Trans. Acousi., Speech, Signal 
Processing, Vol. ASSP-34, no. 6, pp. 1542-1549, December 1986. 

12. Bitmead, R. R. and Anderson B. D. O., "Performance of Adaptive Estimation Al- 
gorithms in Dependent Random Environments," Environments," IEEE Trans. Au- 
tomatic Control, Vol. AC-25, pp. 788-794, August 1980. 

13. Widrow, B. and Stearns, S. D., Adaptive Signal Processing, Prentice-Hall, Inc., 
Englewood Cliffs, New Jersey, 1985. 

14. Wellstead, P. E. and Caldas Pintos J. R., "Self-tuning Filters and Predictors for Two 
Dimensional Systems," International Journal of Control, Vol. 42, no. 2, pp. 457-478, 
1985. 



65 



INITIAL DISTRIBUTION LIST 



No. Copies 

1. Defense Technical Information Center 2 

Cameron Station 

Alexandria, VA 22304-6145 

2. Library, Code 0142 2 

Naval Postgraduate School 

Monterey, CA 93943-5002 

3. Chairman 1 

Electrical and Computer Engineering Department, Code 62 

Naval Postgraduate School 
Monterey, CA 93940-5000 

4. Dr. Murali Tummala 2 

Electrical and Computer Engineering Department, Code 62 

Naval Postgraduate School 
Monterey, CA 93940-5000 

5. Dr. Charles Therrien 1 

Electrical and Computer Engineering Department, Code 62 

Naval Postgraduate School 
Monterey, CA 93940-5000 

6. Professor Will Gersch 1 

Department of Information and Computer Sciences 

Universitv of Hawaii 
Honolulu, HI 96822 

7. Patrol Squadron Forty- Seven 1 

Attn: Lcdr. A. R. Topp 

FPO San Francisco 

San Francisco, CA 96601-5920 

8. Naval Underwater Systems Center 1 

Attn: Lt. R. E. Reinke 

New London, CT 06320 

9. Patrol Squadron Nine 1 

Attn: Lcdr. Steven L. Wilstrup 

FPO San Francisco 

San Francisco, CA 96601-5905 

10. Ocean Systems Center 1 

Commander 

Naval Ocean Svstems Center 
San Diego, CA 92152-5000 



66 



I 

Thesis 

W 644329 Wilstrup 
C -1 Adaptive algorithms 

for two dimensional 
filtering. 



Thecic 

W644329 Wilstrup 

C.l Adaptive algorithms 

for two dimensional 
filtering . 



